program BuoyColl2Bin
!
! Program to put collocated buoy data and scatterometer winds
! into bins and write them to a file that can be plotted
! using the IDL software in genscat/plotting/idl/binning
! The collocated data must be in ASCII format as written by BuoyCollocate
! Output data are also written in ASCII format
!
! Last Modified on: $Date: 2014-08-20 16:35:42 +0200 (Wed, 20 Aug 2014) $
!
!
  use Compiler_Features, only: getarg_genscat, iargc_genscat

  use convert, only: speeddir_to_u, speeddir_to_v

  use binning, only: binned_2D_data_type
  use binning, only: init_2D_bin, add_to_2D_bin, write_2D_bin

  implicit none

  ! data type for collocations
  type :: buoy_coll_type
    integer :: buoy_id
    real    :: buoy_lat
    real    :: buoy_lon
    real    :: buoy_speed
    real    :: buoy_dir
    real    :: speed
    real    :: dir
    logical :: trusted
    integer :: wvc_nr
  end type buoy_coll_type

  type(buoy_coll_type), dimension(:), allocatable :: buoy_coll
  integer, dimension(:), allocatable :: buoy_id ! list of used buoys
  type(binned_2D_data_type) :: binned_data

  integer            :: kbuoy, kbuoy2, nbuoys, kbuoy_id, nbuoy_ids
  integer            :: nr_collocations, nr_buoys, nr_speeds
  integer            :: narg, iarg, error_flag
  real               :: u1, u2, v1, v2, dir1, dir2, lat_tmp, lon_tmp, m_speed
  logical            :: use_model, use_driftcheck, use_stuckcheck, use_smallint, found
  logical            :: free_format
  character(len=256) :: filename_input, filename_output
  character(len=256) :: arg, arg2, tmpstr

  integer            :: buoy_date, buoy_time
  integer            :: scat_date, scat_time
  real               :: scat_lat, scat_lon, scat_speed, scat_dir
  real               :: back_speed, back_dir

  integer,parameter  :: unit_asc = 29

  ! read command line
  narg = iargc_genscat()
  if (narg .eq. 0) then
    print '(A)',' '
    print '(A)','Usage: BuoyColl2Bin -f <ASCII input file> -o <ASCII output file>'
    print '(A)','       [-model] [-nodriftcheck] [-nostuckcheck] [-smallint] [-freeformat]'
    print '(A)',' '
    stop
  endif

  use_model      = .false.
  use_driftcheck = .true.
  use_stuckcheck = .true.
  use_smallint   = .false.
  free_format=.false.

  iarg = 0
  do while (iarg .lt. narg)
    iarg = iarg + 1
    call getarg_genscat(iarg,arg)

    arg2 = arg(2:)
    select case(arg2)

      ! input file
      case('f')
        iarg = iarg + 1
        call getarg_genscat(iarg,filename_input)

      ! output file
      case('o')
        iarg = iarg + 1
        call getarg_genscat(iarg,filename_output)

      ! use model winds
      case('model')
        use_model = .true.

      ! do not check for drifting buoys
      case('nodriftcheck')
        use_driftcheck = .false.

      ! do not check for buoys reporting only low wind speeds
      case('nostuckcheck')
        use_stuckcheck = .false.

      ! use smaller intervals along axes
      case('smallint')
        use_smallint = .true.

      ! free format for reading collocation file
      case ('freeformat')
        free_format=.true.
    end select
  enddo

  ! count number of buoy collocations in file
  error_flag = 0
  nbuoys = 0
  open(unit_asc, file=trim(filename_input), status='old')
  do while (error_flag .eq. 0)
    read(unit_asc, fmt='(A)', iostat=error_flag) tmpstr
    if (error_flag .eq. 0 .and. tmpstr(1:1) .ne. '#') then
      nbuoys = nbuoys + 1
    endif
  enddo
  close(unit_asc)

  allocate(buoy_coll(nbuoys), stat=error_flag)
  if (error_flag .ne. 0) then
    print '(A)',' '
    print '(A,I4)','Error in allocation of buoy_coll, status =', error_flag
    print '(A)',' '
    stop
  endif

  allocate(buoy_id(nbuoys), stat=error_flag)
  if (error_flag .ne. 0) then
    print '(A)',' '
    print '(A,I4)','Error in allocation of buoy_id, status =', error_flag
    print '(A)',' '
    stop
  endif

  ! read buoy data in buoy_coll structure
  ! each buoy measurement can yield 0 or 1 collocations
  nbuoys    = 0
  nbuoy_ids = 0
  open(unit_asc, file=trim(filename_input), status='old')
  do while (error_flag .eq. 0 .and. nbuoys .lt. size(buoy_coll))
    read(unit_asc, fmt='(A)', iostat=error_flag) tmpstr
    if (error_flag .eq. 0 .and. tmpstr(1:1) .ne. '#') then

      ! read formatted file, stop if format seems to be invalid
      nbuoys = nbuoys + 1

      if (free_format) then
      ! unformatted read, stop if error occurs

        if (use_model) then

          ! read buoy winds and model winds
          read(tmpstr, *, iostat=error_flag) &
            buoy_coll(nbuoys)%buoy_id, &
            buoy_coll(nbuoys)%buoy_lat, buoy_coll(nbuoys)%buoy_lon, &
            buoy_date, buoy_time, &
            buoy_coll(nbuoys)%buoy_speed, buoy_coll(nbuoys)%buoy_dir, &
            scat_lat, scat_lon, scat_date, scat_time, &
            scat_speed, scat_dir, &
            buoy_coll(nbuoys)%speed, buoy_coll(nbuoys)%dir

        else

          ! read buoy winds and scatterometer winds
          read(tmpstr, *, iostat=error_flag) &
            buoy_coll(nbuoys)%buoy_id, &
            buoy_coll(nbuoys)%buoy_lat, buoy_coll(nbuoys)%buoy_lon, &
            buoy_date, buoy_time, &
            buoy_coll(nbuoys)%buoy_speed, buoy_coll(nbuoys)%buoy_dir, &
            scat_lat, scat_lon, scat_date, scat_time, &
            buoy_coll(nbuoys)%speed, buoy_coll(nbuoys)%dir, &
            back_speed, back_dir
        endif
        buoy_coll(nbuoys)%trusted = .true.

      else
      ! read formatted file, stop if format seems to be invalid

        if (use_model) then

          ! read buoy winds and model winds
          read(tmpstr, fmt='(I5,F7.2,F8.2,15X,2F6.1,44X,2F6.1)', iostat=error_flag) &
            buoy_coll(nbuoys)%buoy_id, &
            buoy_coll(nbuoys)%buoy_lat, buoy_coll(nbuoys)%buoy_lon, &
            buoy_coll(nbuoys)%buoy_speed, buoy_coll(nbuoys)%buoy_dir, &
            buoy_coll(nbuoys)%speed, buoy_coll(nbuoys)%dir

        else

          ! read buoy winds and scatterometer winds
          read(tmpstr, fmt='(I5,F7.2,F8.2,15X,2F6.1,32X,2F6.1)', iostat=error_flag) &
            buoy_coll(nbuoys)%buoy_id, &
            buoy_coll(nbuoys)%buoy_lat, buoy_coll(nbuoys)%buoy_lon, &
            buoy_coll(nbuoys)%buoy_speed, buoy_coll(nbuoys)%buoy_dir, &
            buoy_coll(nbuoys)%speed, buoy_coll(nbuoys)%dir
        endif
        buoy_coll(nbuoys)%trusted = .true.
      end if

      if (error_flag .ne. 0) then
        print '(A)',' '
        print '(3A,I7)','File ',trim(filename_input), &
          ' seems to contain invalid collocated buoy data near line',nbuoys
        print '(A)',trim(tmpstr)
        print '(A)',' '
        stop
      endif

      ! see if the buoy is already in the list of buoy id's
      found = .false.
      do kbuoy_id = 1, nbuoy_ids
        if (buoy_coll(nbuoys)%buoy_id .eq. buoy_id(kbuoy_id)) then
          found = .true.
          exit
        endif
      enddo
      if (.not. found) then
        nbuoy_ids = nbuoy_ids + 1
        buoy_id(nbuoy_ids) = buoy_coll(nbuoys)%buoy_id
      endif
    endif
  enddo
  close(unit_asc)
  nr_collocations = nbuoys
  nr_buoys        = nbuoy_ids

  ! check if buoy is moored, if it moves more than 0.5 degrees it is considered as drifting
  if (use_driftcheck) then
    buoy_idloop:do kbuoy_id = 1, nbuoy_ids
      lat_tmp = 999.0
      lon_tmp = 999.0
      do kbuoy = 1, nbuoys
        if (buoy_coll(kbuoy)%buoy_id .eq. buoy_id(kbuoy_id)) then
          if (lat_tmp .gt. 100.0) then
            lat_tmp = buoy_coll(kbuoy)%buoy_lat
            lon_tmp = buoy_coll(kbuoy)%buoy_lon
          else
            if (abs(buoy_coll(kbuoy)%buoy_lat - lat_tmp) .gt. 0.5 .or. &
                abs(buoy_coll(kbuoy)%buoy_lon - lon_tmp) .gt. 0.5) then
              print '(A,I6,A)','Buoy', buoy_id(kbuoy_id), ' seems to be drifting, it will not be used'
              nr_buoys = nr_buoys - 1

              ! flag all buoy collocations of drifting buoy
              do kbuoy2 = 1, nbuoys
                if (buoy_coll(kbuoy2)%buoy_id .eq. buoy_id(kbuoy_id)) then
                  nr_collocations = nr_collocations - 1
                  buoy_coll(kbuoy2)%trusted = .false.
                endif
              enddo
              cycle buoy_idloop
            endif
          endif
        endif
      enddo
    enddo buoy_idloop
  endif

  ! check if buoy reports only low wind speeds, its anemometer may be stuck
  if (use_stuckcheck) then
    buoy_idloop2:do kbuoy_id = 1, nbuoy_ids
      nr_speeds = 0
      m_speed   = 0.0
      do kbuoy = 1, nbuoys
        if (buoy_coll(kbuoy)%buoy_id .eq. buoy_id(kbuoy_id)) then
          nr_speeds = nr_speeds + 1
          m_speed   = m_speed + buoy_coll(kbuoy)%buoy_speed
        endif
      enddo
      if (nr_speeds .gt. 0) m_speed = m_speed / real(nr_speeds)

      if (nr_speeds .gt. 2 .and. m_speed .lt. 1.5) then  ! mean wind speed below 1.5 m/s
        print '(A,I6,A)','Buoy', buoy_id(kbuoy_id), ' only reports low wind speeds, it will not be used'
        nr_buoys = nr_buoys - 1

        ! flag all buoy collocations of drifting buoy
        do kbuoy2 = 1, nbuoys
          if (buoy_coll(kbuoy2)%buoy_id .eq. buoy_id(kbuoy_id)) then
            nr_collocations = nr_collocations - 1
            buoy_coll(kbuoy2)%trusted = .false.
          endif
        enddo
      endif

    enddo buoy_idloop2
  endif

  ! skip certain WVC numbers if needed
! do kbuoy = 1, nbuoys
!   if (buoy_coll(kbuoy)%wvc_nr .le. 4 .or. buoy_coll(kbuoy)%wvc_nr .ge. 33) then
!     buoy_coll(kbuoy)%trusted = .false.
!   endif
! enddo

  ! open ASCII output file and write heading
  open(unit_asc, file=filename_output)
  if (use_model) then
    write(unit_asc,'(A)') 'model'
  else
    write(unit_asc,'(A)') 'scatterometer'
  endif
  write(unit_asc,'(A)') 'buoy'

  ! calculate and write binned data of u component
  ! add 0.01 to direction to prevent binning problems around 180 and 360 degrees
  if (use_smallint) then
    call init_2D_bin(binned_data, 100, -25.0, 25.0)
  else
    call init_2D_bin(binned_data, 50, -25.0, 25.0)
  endif
  do kbuoy = 1, nbuoys
    if (buoy_coll(kbuoy)%trusted) then
      u1 = speeddir_to_u(buoy_coll(kbuoy)%buoy_speed, buoy_coll(kbuoy)%buoy_dir+0.01)
      u2 = speeddir_to_u(buoy_coll(kbuoy)%speed, buoy_coll(kbuoy)%dir+0.01)
      call add_to_2D_bin(binned_data, u1, u2)
    endif
  enddo
  call write_2D_bin(binned_data, unit_asc, silent=.true.)

  ! calculate and write binned data of v component
  ! add 0.01 to direction to prevent binning problems around 180 and 360 degrees
  if (use_smallint) then
    call init_2D_bin(binned_data, 100, -25.0, 25.0)
  else
    call init_2D_bin(binned_data, 50, -25.0, 25.0)
  endif
  do kbuoy = 1, nbuoys
    if (buoy_coll(kbuoy)%trusted) then
      v1 = speeddir_to_v(buoy_coll(kbuoy)%buoy_speed, buoy_coll(kbuoy)%buoy_dir+0.01)
      v2 = speeddir_to_v(buoy_coll(kbuoy)%speed, buoy_coll(kbuoy)%dir+0.01)
      call add_to_2D_bin(binned_data, v1, v2)
    endif
  enddo
  call write_2D_bin(binned_data, unit_asc, silent=.true.)

  ! calculate and write binned data of wind speed
  if (use_smallint) then
    call init_2D_bin(binned_data, 50, 0.0, 25.0)
  else
    call init_2D_bin(binned_data, 25, 0.0, 25.0)
  endif
  do kbuoy = 1, nbuoys
    if (buoy_coll(kbuoy)%trusted) then
      call add_to_2D_bin(binned_data, buoy_coll(kbuoy)%buoy_speed, buoy_coll(kbuoy)%speed)
    endif
  enddo
  call write_2D_bin(binned_data, unit_asc, silent=.true.)

  ! calculate and write binned data of wind direction, only for winds >= 4 m/s
  ! the buoy wind directions in MARS are rounded to 10 degrees
  ! round buoy and scat directions to 10 degrees to prevent binning problems
  ! subtract 0.01 from scat directions to prevent 2.5 degrees bias for MSS data
  !
  ! directions are rounded to 2.5 degrees in other cases to prevent problems
  ! with MSS data which have a direction range of 0.0-357.5 rather than 0.0-360.0
  if (use_smallint) then
    call init_2D_bin(binned_data, 144, 0.0, 360.0)
  else
    call init_2D_bin(binned_data, 36, 0.0, 360.0)
  endif
  do kbuoy = 1, nbuoys
    if (buoy_coll(kbuoy)%trusted) then
      if (buoy_coll(kbuoy)%buoy_speed .ge. 4.0) then
        if (use_smallint) then
          dir1 = real(nint(buoy_coll(kbuoy)%buoy_dir / 2.5)) * 2.5
          dir2 = real(nint(buoy_coll(kbuoy)%dir / 2.5)) * 2.5
        else
          dir1 = real(nint(buoy_coll(kbuoy)%buoy_dir / 10.0)) * 10.0
          dir2 = real(nint((buoy_coll(kbuoy)%dir-0.01) / 10.0)) * 10.0
        endif
        call add_to_2D_bin(binned_data, dir1, dir2)
      endif
    endif
  enddo
  call write_2D_bin(binned_data, unit_asc, silent=.true.)

  ! close output file
  close(unit_asc)

  if (allocated(buoy_coll)) deallocate(buoy_coll)
  if (allocated(buoy_id)) deallocate(buoy_id)

  ! report number of collocations and number of buoys
  print '(A,I7,A,I6,A)','Found', nr_collocations, ' collocations from', nr_buoys, &
                        ' different buoys'

end program BuoyColl2Bin
